library(data.table)
library(Seurat)
library(harmony)
library(SingleR)
library(tidyseurat)
library(tidyverse)
source('00_util_scripts/mod_seurat.R')

# read raw data -------------
count_mtx <- fread('mission/FPP/psoriasis/GSE162183_Raw_gene_counts_matrix.tab.gz')

count_mtx[1:5, 1:5]

column_to_rownames(count_mtx, 'V1') -> count_mtx

sobj <- CreateSeuratObject(count_mtx, min.cells = 3, min.features = 200)

sobj <- sobj |>
  PercentageFeatureSet('^MT-', col.name = 'mito.ratio')

sobj |> VlnPlot('mito.ratio', pt.size = 0)

sobj <- sobj |>
  filter(mito.ratio < 10)

sobj <- sobj |>
  mutate(group = str_extract(orig.ident, 'Ctrl|Psor'))

# pre-process -------------
sobj <- sobj |>
  quick_process_seurat()

# identify cell type ---------
hpca <- celldex::HumanPrimaryCellAtlasData()

sobj <- sobj |>
  mark_cell_type_singler(hpca, new_label = 'hpca_main')

# save and read -----------
write_rds(sobj, 'mission/FPP/psoriasis/gao2021harbin.rds')

sobj <- read_rds('mission/FPP/psoriasis/gao2021harbin.rds')

sobj |>
  filter(str_detect(hpca_main, 'Kera')) |>
  write_rds('mission/FPP/psoriasis/gao2021kera.rds')

sobj.kera <- read_rds('mission/FPP/psoriasis/gao2021kera.rds')

# visualization mission -----------
kegg_mva <-
  c('Acat1','Acat2','Hmgcs1','Hmgcs2','Hmgcr','Mvk','Pmvk','Mvd','Idi1','Idi2','Fdps')|>
  str_to_upper()

mva_fct <- tibble(gene = kegg_mva, ordered = fct_inorder(kegg_mva))

kegg_d.mva <-
  c('FDFT1','LSS','SQLE','GGPS1','DOLK','PDSS1','PDSS2') |>
  fct_inorder()

key_cytokine <- c('Il6','Ccl20','Tslp','Flt3l','Csf2','Tnf') |>
  str_to_upper()

celltyps.list <- sobj$manual.main |> unique()

celltyps.list <- celltyps.list |>
  set_names(celltyps.list)

Idents(sobj) <- 'manual.main'

celltyps.psor.vs.ctrl <- celltyps.list |>
  map(\(x)sobj |>
        FindMarkers(ident.1 = 'Psor', group.by = 'group', subset.ident = x) |>
        as_tibble(rownames = 'gene'), .progress = T) |>
  list_rbind(names_to = 'cell.type')

celltyps.psor.vs.ctrl |>
  right_join(mva_fct) |>
  filter(!is.na(avg_log2FC)) |>
  ggplot(aes(ordered, cell.type, color = avg_log2FC, size = -log10(p_val_adj))) +
  geom_point() +
  scale_size(range = c(0,3)) +
  scale_color_distiller(palette = 'RdYlBu') +
  labs(x = 'gene', title = 'MVA pathway in psoriasis vs control') +
  theme_jpub +
  rotate_x_text(45)

publish_pdf('mission/FPP/psoriasis/psovhc.mva.allskin.pdf', width = 70)

celltyps.psor.vs.ctrl |>
  inner_join(mva_fct) |>
  write_csv('mission/FPP/pub_source_data/fig.N.MVA.psoriasis.vs.HC.csv')

celltyps.psor.vs.ctrl |>
  filter(gene %in% key_cytokine) |>
  ggplot(aes(gene, cell.type, color = avg_log2FC, size = -log10(p_val_adj))) +
  geom_point() +
  scale_color_distiller(palette = 'RdBu',limits = c(-4,4)) +
  theme_pubr(x.text.angle = 45) +
  labs(x = 'gene')

sobj.kera |>
  get_abundance_sc_long('CCL20') |>
  mutate(group = ifelse(str_detect(.cell, 'Ctrl'), 'Ctrl', 'Psoriasis')) |>
  ggplot(aes(group, .abundance_RNA, fill = group)) +
  geom_violin() +
  stat_summary(geom = 'crossbar', fun = 'mean',
               color = 'red', width = .3) +
  theme_pubr(legend = 'none') +
  labs(title = 'CCL20 RNA expression in keratinocytes',
       y = 'Normalized expression')

sobj.kera <- sobj.kera |>
  quick_process_seurat(skip_norm = T)

sobj.kera |>
  FindAllMarkers(features = 'TRPV3', only.pos = T)

sobj.kera <- sobj.kera |>
  mutate(trpv3.hi = ifelse(seurat_clusters %in% c(2,7,6,8),
                           'TRPV3_high', 'TRPV3_low'))

sobj.kera |>
  VlnPlot(features = 'TRPV3', group.by = 'trpv3.hi', pt.size = 0) +
  stat_summary(geom = 'crossbar', fun = 'mean',
               color = 'red', width = .3) +
  theme_pubr(legend = 'none') +
  labs(title = 'CCL20 RNA expression in keratinocytes',
       y = 'Normalized expression', x = 'Subtype')

sobj.kera |>
  DimPlot(group.by = 'trpv3.hi') +
  ggtitle('Keratinocyte subtype on TRPV3 expression')

sobj.kera <- sobj.kera |>
  mutate(subgroup = str_c(trpv3.hi,'_', group)) 

sobj.kera |>
  DotPlot(kegg_mva, group.by = 'subgroup') +
  theme_pubr(x.text.angle = 45, legend = 'right') +
  labs(y = 'Subgroup', x = 'Gene',
       title = 'MVA pathway in TRPV3-hi/lo keratinocyte in psoriasis/ctrl skin')

sobj.kera |>
  DotPlot(key_cytokine, group.by = 'subgroup') +
  theme_pubr(x.text.angle = 45, legend = 'right') +
  labs(y = 'Subgroup', x = 'Gene',
       title = 'Cytokine in TRPV3-hi/lo keratinocyte in psoriasis/ctrl skin')

## V3h-KC in skin -------------
v3h.kc <- sobj.kera |>
  filter(trpv3.hi == 'TRPV3_high') |>
  colnames()

v3h.kc |> head()

sobj <- sobj |>
  mutate(manual.main = case_when(.cell %in% v3h.kc ~ 'V3-hi KC',
                                 str_detect(hpca_main, 'Kera') ~ 'V3-lo KC',
                                 .default = hpca_main))

sobj |> DimPlot(group.by = 'manual.main')

Idents(sobj) <- 'manual.main'

sobj |> write_rds('mission/FPP/psoriasis/gao2021harbin.rds')

pso.fine <- sobj$manual.main |> unique()

skin.fc.psovhc <- pso.fine |>
  map(\(x)sobj |> FindMarkers(group.by = 'group', ident.1 = 'Psor',
                              subset.ident = x) |>
        as_tibble(rownames = 'gene') |>
        mutate(celltype = x),
      .progress = T) |>
  list_rbind()

skin.fc.psovhc |>
  filter(p_val_adj < .05) |>
  write_csv('mission/FPP/psoriasis/skin.psovhc.alltype.deg.csv')

skin.fc.psovhc <-
  read_csv('mission/FPP/psoriasis/skin.psovhc.alltype.deg.csv')

### UPR chaperone -----------
upr.chap.ref <-
  read_csv('mission/FPP/UPR.associated.chaperone.gene.csv')

skin.fc.psovhc |>
  filter(gene %in% upr.chap.ref$SYMBOL) |>
  ggplot(aes(y = celltype, x = gene, size = -log10(p_val_adj),
             color = avg_log2FC)) +
  geom_point() +
  scale_color_distiller(palette = 'RdYlBu') +
  scale_size(range = c(0,3)) +
  theme_pubr(legend = 'right') +
  theme_jpub +
  labs(x = 'Gene', y = 'Cell type',
       title = 'ER stress-associated chaperon: Psoriasis vs Ctrl') +
  RotatedAxis()

publish_pdf('mission/FPP/psoriasis/psovhc.skin.chaperone.logfc.pdf', width = 80)

### GO GSEA -------------
skin.psovhc.gogse <- pso.fine |>
  map(\(x){skin.fc.psovhc |>
      filter(celltype == x) |>
      pull(avg_log2FC, name = gene) |>
      sort(decreasing = T) |>
      gseGO(ont = 'ALL', OrgDb = 'org.Hs.eg.db',
            keyType = 'SYMBOL',pvalueCutoff = 1) |>
      pluck('result') |>
      mutate(cell.type = x)}) |>
  list_rbind() |>
  as_tibble()

skin.psovhc.gogse |>
  select(-c(setSize, enrichmentScore)) |>
  write_csv('mission/FPP/psoriasis/skin.all.psovhc.go.gsea.csv')

skin.psovhc.gogse <-
  read_csv('mission/FPP/psoriasis/skin.all.psovhc.go.gsea.csv')

plot.gsea.dot <- function(df) {
  df |> ggplot(aes(.data$cell.type, str_wrap(.data$Description, 50),
                   size = -log10(.data$pvalue), color = .data$NES)) +
    geom_point() +
    scale_size(range = c(0,3)) +
    scale_color_distiller(palette = 'RdYlBu') +
    theme_jpub +
    rotate_x_text(45) +
    labs(title = 'Differential pathway enrichment in psoriasis vs ctrl',
         y = 'GO term', x = 'Cell type')
}

skin.psovhc.gogse |>
  filter(str_detect(Description, 'reactive oxygen')) |>
  plot.gsea.dot()

publish_pdf('mission/FPP/psoriasis/skin.psovhc.ros.gsea.pdf', width = 77)

skin.psovhc.gogse |>
  filter(str_detect(Description, 'TOR')) |>
  plot.gsea.dot()

publish_pdf('mission/FPP/psoriasis/skin.psovhc.TOR.gsea.pdf', width = 70)

skin.psovhc.gogse |>
  filter(str_detect(Description, '[^o]autophag')) |>
  plot.gsea.dot()

publish_pdf('mission/FPP/psoriasis/skin.psovhc.autophagy.gsea.pdf', width = 77)

skin.psovhc.gogse |>
  filter(str_detect(Description, 'pyropto|inflammasome')) |>
  plot.gsea.dot()

publish_pdf('mission/FPP/psoriasis/skin.psovhc.inflammasome.gsea.pdf',
            width = 75)

skin.psovhc.gogse |>
  filter(str_detect(Description, 'toll')) |>
  plot.gsea.dot()

publish_pdf('mission/FPP/psoriasis/skin.psovhc.TLR.gsea.pdf',
            width = 75)

skin.psovhc.gogse |>
  filter(str_detect(Description, 'interleukin-(1$|1 )')) |>
  plot.gsea.dot()

publish_pdf('mission/FPP/psoriasis/skin.psovhc.IL1.gsea.pdf',
            width = 77)
